1 // Written in the D programming language. 2 3 /** 4 * Contains the elementary mathematical functions (powers, roots, 5 * and trigonometric functions), and low-level floating-point operations. 6 * Mathematical special functions are available in `std.mathspecial`. 7 * 8 $(SCRIPT inhibitQuickIndex = 1;) 9 10 $(DIVC quickindex, 11 $(BOOKTABLE , 12 $(TR $(TH Category) $(TH Members) ) 13 $(TR $(TDNW $(SUBMODULE Constants, constants)) $(TD 14 $(SUBREF constants, E) 15 $(SUBREF constants, PI) 16 $(SUBREF constants, PI_2) 17 $(SUBREF constants, PI_4) 18 $(SUBREF constants, M_1_PI) 19 $(SUBREF constants, M_2_PI) 20 $(SUBREF constants, M_2_SQRTPI) 21 $(SUBREF constants, LN10) 22 $(SUBREF constants, LN2) 23 $(SUBREF constants, LOG2) 24 $(SUBREF constants, LOG2E) 25 $(SUBREF constants, LOG2T) 26 $(SUBREF constants, LOG10E) 27 $(SUBREF constants, SQRT2) 28 $(SUBREF constants, SQRT1_2) 29 )) 30 $(TR $(TDNW $(SUBMODULE Algebraic, algebraic)) $(TD 31 $(SUBREF algebraic, abs) 32 $(SUBREF algebraic, fabs) 33 $(SUBREF algebraic, sqrt) 34 $(SUBREF algebraic, cbrt) 35 $(SUBREF algebraic, hypot) 36 $(SUBREF algebraic, poly) 37 $(SUBREF algebraic, nextPow2) 38 $(SUBREF algebraic, truncPow2) 39 )) 40 $(TR $(TDNW $(SUBMODULE Trigonometry, trigonometry)) $(TD 41 $(SUBREF trigonometry, sin) 42 $(SUBREF trigonometry, cos) 43 $(SUBREF trigonometry, tan) 44 $(SUBREF trigonometry, asin) 45 $(SUBREF trigonometry, acos) 46 $(SUBREF trigonometry, atan) 47 $(SUBREF trigonometry, atan2) 48 $(SUBREF trigonometry, sinh) 49 $(SUBREF trigonometry, cosh) 50 $(SUBREF trigonometry, tanh) 51 $(SUBREF trigonometry, asinh) 52 $(SUBREF trigonometry, acosh) 53 $(SUBREF trigonometry, atanh) 54 )) 55 $(TR $(TDNW $(SUBMODULE Rounding, rounding)) $(TD 56 $(SUBREF rounding, ceil) 57 $(SUBREF rounding, floor) 58 $(SUBREF rounding, round) 59 $(SUBREF rounding, lround) 60 $(SUBREF rounding, trunc) 61 $(SUBREF rounding, rint) 62 $(SUBREF rounding, lrint) 63 $(SUBREF rounding, nearbyint) 64 $(SUBREF rounding, rndtol) 65 $(SUBREF rounding, quantize) 66 )) 67 $(TR $(TDNW $(SUBMODULE Exponentiation & Logarithms, exponential)) $(TD 68 $(SUBREF exponential, pow) 69 $(SUBREF exponential, powmod) 70 $(SUBREF exponential, exp) 71 $(SUBREF exponential, exp2) 72 $(SUBREF exponential, expm1) 73 $(SUBREF exponential, ldexp) 74 $(SUBREF exponential, frexp) 75 $(SUBREF exponential, log) 76 $(SUBREF exponential, log2) 77 $(SUBREF exponential, log10) 78 $(SUBREF exponential, logb) 79 $(SUBREF exponential, ilogb) 80 $(SUBREF exponential, log1p) 81 $(SUBREF exponential, scalbn) 82 )) 83 $(TR $(TDNW $(SUBMODULE Remainder, remainder)) $(TD 84 $(SUBREF remainder, fmod) 85 $(SUBREF remainder, modf) 86 $(SUBREF remainder, remainder) 87 $(SUBREF remainder, remquo) 88 )) 89 $(TR $(TDNW $(SUBMODULE Floating-point operations, operations)) $(TD 90 $(SUBREF operations, approxEqual) 91 $(SUBREF operations, feqrel) 92 $(SUBREF operations, fdim) 93 $(SUBREF operations, fmax) 94 $(SUBREF operations, fmin) 95 $(SUBREF operations, fma) 96 $(SUBREF operations, isClose) 97 $(SUBREF operations, nextDown) 98 $(SUBREF operations, nextUp) 99 $(SUBREF operations, nextafter) 100 $(SUBREF operations, NaN) 101 $(SUBREF operations, getNaNPayload) 102 $(SUBREF operations, cmp) 103 )) 104 $(TR $(TDNW $(SUBMODULE Introspection, traits)) $(TD 105 $(SUBREF traits, isFinite) 106 $(SUBREF traits, isIdentical) 107 $(SUBREF traits, isInfinity) 108 $(SUBREF traits, isNaN) 109 $(SUBREF traits, isNormal) 110 $(SUBREF traits, isSubnormal) 111 $(SUBREF traits, signbit) 112 $(SUBREF traits, sgn) 113 $(SUBREF traits, copysign) 114 $(SUBREF traits, isPowerOf2) 115 )) 116 $(TR $(TDNW $(SUBMODULE Hardware Control, hardware)) $(TD 117 $(SUBREF hardware, IeeeFlags) 118 $(SUBREF hardware, ieeeFlags) 119 $(SUBREF hardware, resetIeeeFlags) 120 $(SUBREF hardware, FloatingPointControl) 121 )) 122 ) 123 ) 124 125 * The functionality closely follows the IEEE754-2008 standard for 126 * floating-point arithmetic, including the use of camelCase names rather 127 * than C99-style lower case names. All of these functions behave correctly 128 * when presented with an infinity or NaN. 129 * 130 * The following IEEE 'real' formats are currently supported: 131 * $(UL 132 * $(LI 64 bit Big-endian 'double' (eg PowerPC)) 133 * $(LI 128 bit Big-endian 'quadruple' (eg SPARC)) 134 * $(LI 64 bit Little-endian 'double' (eg x86-SSE2)) 135 * $(LI 80 bit Little-endian, with implied bit 'real80' (eg x87, Itanium)) 136 * $(LI 128 bit Little-endian 'quadruple' (not implemented on any known processor!)) 137 * $(LI Non-IEEE 128 bit Big-endian 'doubledouble' (eg PowerPC) has partial support) 138 * ) 139 * Unlike C, there is no global 'errno' variable. Consequently, almost all of 140 * these functions are pure nothrow. 141 * 142 * Macros: 143 * SUBMODULE = $(MREF_ALTTEXT $1, std, math, $2) 144 * SUBREF = $(REF_ALTTEXT $(TT $2), $2, std, math, $1)$(NBSP) 145 * 146 * Copyright: Copyright The D Language Foundation 2000 - 2011. 147 * D implementations of tan, atan, atan2, exp, expm1, exp2, log, log10, log1p, 148 * log2, floor, ceil and lrint functions are based on the CEPHES math library, 149 * which is Copyright (C) 2001 Stephen L. Moshier $(LT)steve@moshier.net$(GT) 150 * and are incorporated herein by permission of the author. The author 151 * reserves the right to distribute this material elsewhere under different 152 * copying permissions. These modifications are distributed here under 153 * the following terms: 154 * License: $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 155 * Authors: $(HTTP digitalmars.com, Walter Bright), Don Clugston, 156 * Conversion of CEPHES math library to D by Iain Buclaw and David Nadlinger 157 * Source: $(PHOBOSSRC std/math/package.d) 158 */ 159 module std.math; 160 161 public import std.math.algebraic; 162 public import std.math.constants; 163 public import std.math.exponential; 164 public import std.math.operations; 165 public import std.math.hardware; 166 public import std.math.remainder; 167 public import std.math.rounding; 168 public import std.math.traits; 169 public import std.math.trigonometry; 170 171 // @@@DEPRECATED_2.102@@@ 172 // Note: Exposed accidentally, should be deprecated / removed 173 deprecated("std.meta.AliasSeq was unintentionally available from std.math " 174 ~ "and will be removed after 2.102. Please import std.meta instead") 175 public import std.meta : AliasSeq; 176 177 package(std): // Not public yet 178 /* Return the value that lies halfway between x and y on the IEEE number line. 179 * 180 * Formally, the result is the arithmetic mean of the binary significands of x 181 * and y, multiplied by the geometric mean of the binary exponents of x and y. 182 * x and y must have the same sign, and must not be NaN. 183 * Note: this function is useful for ensuring O(log n) behaviour in algorithms 184 * involving a 'binary chop'. 185 * 186 * Special cases: 187 * If x and y are within a factor of 2, (ie, feqrel(x, y) > 0), the return value 188 * is the arithmetic mean (x + y) / 2. 189 * If x and y are even powers of 2, the return value is the geometric mean, 190 * ieeeMean(x, y) = sqrt(x * y). 191 * 192 */ 193 T ieeeMean(T)(const T x, const T y) @trusted pure nothrow @nogc 194 in 195 { 196 // both x and y must have the same sign, and must not be NaN. 197 assert(signbit(x) == signbit(y)); 198 assert(x == x && y == y); 199 } 200 do 201 { 202 // Runtime behaviour for contract violation: 203 // If signs are opposite, or one is a NaN, return 0. 204 if (!((x >= 0 && y >= 0) || (x <= 0 && y <= 0))) return 0.0; 205 206 // The implementation is simple: cast x and y to integers, 207 // average them (avoiding overflow), and cast the result back to a floating-point number. 208 209 alias F = floatTraits!(T); 210 T u; 211 static if (F.realFormat == RealFormat.ieeeExtended || 212 F.realFormat == RealFormat.ieeeExtended53) 213 { 214 // There's slight additional complexity because they are actually 215 // 79-bit reals... 216 ushort *ue = cast(ushort *)&u; 217 ulong *ul = cast(ulong *)&u; 218 ushort *xe = cast(ushort *)&x; 219 ulong *xl = cast(ulong *)&x; 220 ushort *ye = cast(ushort *)&y; 221 ulong *yl = cast(ulong *)&y; 222 223 // Ignore the useless implicit bit. (Bonus: this prevents overflows) 224 ulong m = ((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL); 225 226 // @@@ BUG? @@@ 227 // Cast shouldn't be here 228 ushort e = cast(ushort) ((xe[F.EXPPOS_SHORT] & F.EXPMASK) 229 + (ye[F.EXPPOS_SHORT] & F.EXPMASK)); 230 if (m & 0x8000_0000_0000_0000L) 231 { 232 ++e; 233 m &= 0x7FFF_FFFF_FFFF_FFFFL; 234 } 235 // Now do a multi-byte right shift 236 const uint c = e & 1; // carry 237 e >>= 1; 238 m >>>= 1; 239 if (c) 240 m |= 0x4000_0000_0000_0000L; // shift carry into significand 241 if (e) 242 *ul = m | 0x8000_0000_0000_0000L; // set implicit bit... 243 else 244 *ul = m; // ... unless exponent is 0 (subnormal or zero). 245 246 ue[4]= e | (xe[F.EXPPOS_SHORT]& 0x8000); // restore sign bit 247 } 248 else static if (F.realFormat == RealFormat.ieeeQuadruple) 249 { 250 // This would be trivial if 'ucent' were implemented... 251 ulong *ul = cast(ulong *)&u; 252 ulong *xl = cast(ulong *)&x; 253 ulong *yl = cast(ulong *)&y; 254 255 // Multi-byte add, then multi-byte right shift. 256 import core.checkedint : addu; 257 bool carry; 258 ulong ml = addu(xl[MANTISSA_LSB], yl[MANTISSA_LSB], carry); 259 260 ulong mh = carry + (xl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL) + 261 (yl[MANTISSA_MSB] & 0x7FFF_FFFF_FFFF_FFFFL); 262 263 ul[MANTISSA_MSB] = (mh >>> 1) | (xl[MANTISSA_MSB] & 0x8000_0000_0000_0000); 264 ul[MANTISSA_LSB] = (ml >>> 1) | (mh & 1) << 63; 265 } 266 else static if (F.realFormat == RealFormat.ieeeDouble) 267 { 268 ulong *ul = cast(ulong *)&u; 269 ulong *xl = cast(ulong *)&x; 270 ulong *yl = cast(ulong *)&y; 271 ulong m = (((*xl) & 0x7FFF_FFFF_FFFF_FFFFL) 272 + ((*yl) & 0x7FFF_FFFF_FFFF_FFFFL)) >>> 1; 273 m |= ((*xl) & 0x8000_0000_0000_0000L); 274 *ul = m; 275 } 276 else static if (F.realFormat == RealFormat.ieeeSingle) 277 { 278 uint *ul = cast(uint *)&u; 279 uint *xl = cast(uint *)&x; 280 uint *yl = cast(uint *)&y; 281 uint m = (((*xl) & 0x7FFF_FFFF) + ((*yl) & 0x7FFF_FFFF)) >>> 1; 282 m |= ((*xl) & 0x8000_0000); 283 *ul = m; 284 } 285 else 286 { 287 assert(0, "Not implemented"); 288 } 289 return u; 290 } 291 292 @safe pure nothrow @nogc unittest 293 { 294 assert(ieeeMean(-0.0,-1e-20)<0); 295 assert(ieeeMean(0.0,1e-20)>0); 296 297 assert(ieeeMean(1.0L,4.0L)==2L); 298 assert(ieeeMean(2.0*1.013,8.0*1.013)==4*1.013); 299 assert(ieeeMean(-1.0L,-4.0L)==-2L); 300 assert(ieeeMean(-1.0,-4.0)==-2); 301 assert(ieeeMean(-1.0f,-4.0f)==-2f); 302 assert(ieeeMean(-1.0,-2.0)==-1.5); 303 assert(ieeeMean(-1*(1+8*real.epsilon),-2*(1+8*real.epsilon)) 304 ==-1.5*(1+5*real.epsilon)); 305 assert(ieeeMean(0x1p60,0x1p-10)==0x1p25); 306 307 static if (floatTraits!(real).realFormat == RealFormat.ieeeExtended) 308 { 309 assert(ieeeMean(1.0L,real.infinity)==0x1p8192L); 310 assert(ieeeMean(0.0L,real.infinity)==1.5); 311 } 312 assert(ieeeMean(0.5*real.min_normal*(1-4*real.epsilon),0.5*real.min_normal) 313 == 0.5*real.min_normal*(1-2*real.epsilon)); 314 } 315 316 317 // The following IEEE 'real' formats are currently supported. 318 version (LittleEndian) 319 { 320 static assert(real.mant_dig == 53 || real.mant_dig == 64 321 || real.mant_dig == 113, 322 "Only 64-bit, 80-bit, and 128-bit reals"~ 323 " are supported for LittleEndian CPUs"); 324 } 325 else 326 { 327 static assert(real.mant_dig == 53 || real.mant_dig == 113, 328 "Only 64-bit and 128-bit reals are supported for BigEndian CPUs."); 329 } 330 331 // Underlying format exposed through floatTraits 332 enum RealFormat 333 { 334 ieeeHalf, 335 ieeeSingle, 336 ieeeDouble, 337 ieeeExtended, // x87 80-bit real 338 ieeeExtended53, // x87 real rounded to precision of double. 339 ibmExtended, // IBM 128-bit extended 340 ieeeQuadruple, 341 } 342 343 // Constants used for extracting the components of the representation. 344 // They supplement the built-in floating point properties. 345 template floatTraits(T) 346 { 347 import std.traits : Unqual; 348 349 // EXPMASK is a ushort mask to select the exponent portion (without sign) 350 // EXPSHIFT is the number of bits the exponent is left-shifted by in its ushort 351 // EXPBIAS is the exponent bias - 1 (exp == EXPBIAS yields ×2^-1). 352 // EXPPOS_SHORT is the index of the exponent when represented as a ushort array. 353 // SIGNPOS_BYTE is the index of the sign when represented as a ubyte array. 354 // RECIP_EPSILON is the value such that (smallest_subnormal) * RECIP_EPSILON == T.min_normal 355 enum Unqual!T RECIP_EPSILON = (1/T.epsilon); 356 static if (T.mant_dig == 24) 357 { 358 // Single precision float 359 enum ushort EXPMASK = 0x7F80; 360 enum ushort EXPSHIFT = 7; 361 enum ushort EXPBIAS = 0x3F00; 362 enum uint EXPMASK_INT = 0x7F80_0000; 363 enum uint MANTISSAMASK_INT = 0x007F_FFFF; 364 enum realFormat = RealFormat.ieeeSingle; 365 version (LittleEndian) 366 { 367 enum EXPPOS_SHORT = 1; 368 enum SIGNPOS_BYTE = 3; 369 } 370 else 371 { 372 enum EXPPOS_SHORT = 0; 373 enum SIGNPOS_BYTE = 0; 374 } 375 } 376 else static if (T.mant_dig == 53) 377 { 378 static if (T.sizeof == 8) 379 { 380 // Double precision float, or real == double 381 enum ushort EXPMASK = 0x7FF0; 382 enum ushort EXPSHIFT = 4; 383 enum ushort EXPBIAS = 0x3FE0; 384 enum uint EXPMASK_INT = 0x7FF0_0000; 385 enum uint MANTISSAMASK_INT = 0x000F_FFFF; // for the MSB only 386 enum ulong MANTISSAMASK_LONG = 0x000F_FFFF_FFFF_FFFF; 387 enum realFormat = RealFormat.ieeeDouble; 388 version (LittleEndian) 389 { 390 enum EXPPOS_SHORT = 3; 391 enum SIGNPOS_BYTE = 7; 392 } 393 else 394 { 395 enum EXPPOS_SHORT = 0; 396 enum SIGNPOS_BYTE = 0; 397 } 398 } 399 else static if (T.sizeof == 12) 400 { 401 // Intel extended real80 rounded to double 402 enum ushort EXPMASK = 0x7FFF; 403 enum ushort EXPSHIFT = 0; 404 enum ushort EXPBIAS = 0x3FFE; 405 enum realFormat = RealFormat.ieeeExtended53; 406 version (LittleEndian) 407 { 408 enum EXPPOS_SHORT = 4; 409 enum SIGNPOS_BYTE = 9; 410 } 411 else 412 { 413 enum EXPPOS_SHORT = 0; 414 enum SIGNPOS_BYTE = 0; 415 } 416 } 417 else 418 static assert(false, "No traits support for " ~ T.stringof); 419 } 420 else static if (T.mant_dig == 64) 421 { 422 // Intel extended real80 423 enum ushort EXPMASK = 0x7FFF; 424 enum ushort EXPSHIFT = 0; 425 enum ushort EXPBIAS = 0x3FFE; 426 enum realFormat = RealFormat.ieeeExtended; 427 version (LittleEndian) 428 { 429 enum EXPPOS_SHORT = 4; 430 enum SIGNPOS_BYTE = 9; 431 } 432 else 433 { 434 enum EXPPOS_SHORT = 0; 435 enum SIGNPOS_BYTE = 0; 436 } 437 } 438 else static if (T.mant_dig == 113) 439 { 440 // Quadruple precision float 441 enum ushort EXPMASK = 0x7FFF; 442 enum ushort EXPSHIFT = 0; 443 enum ushort EXPBIAS = 0x3FFE; 444 enum realFormat = RealFormat.ieeeQuadruple; 445 version (LittleEndian) 446 { 447 enum EXPPOS_SHORT = 7; 448 enum SIGNPOS_BYTE = 15; 449 } 450 else 451 { 452 enum EXPPOS_SHORT = 0; 453 enum SIGNPOS_BYTE = 0; 454 } 455 } 456 else static if (T.mant_dig == 106) 457 { 458 // IBM Extended doubledouble 459 enum ushort EXPMASK = 0x7FF0; 460 enum ushort EXPSHIFT = 4; 461 enum realFormat = RealFormat.ibmExtended; 462 463 // For IBM doubledouble the larger magnitude double comes first. 464 // It's really a double[2] and arrays don't index differently 465 // between little and big-endian targets. 466 enum DOUBLEPAIR_MSB = 0; 467 enum DOUBLEPAIR_LSB = 1; 468 469 // The exponent/sign byte is for most significant part. 470 version (LittleEndian) 471 { 472 enum EXPPOS_SHORT = 3; 473 enum SIGNPOS_BYTE = 7; 474 } 475 else 476 { 477 enum EXPPOS_SHORT = 0; 478 enum SIGNPOS_BYTE = 0; 479 } 480 } 481 else 482 static assert(false, "No traits support for " ~ T.stringof); 483 } 484 485 // These apply to all floating-point types 486 version (LittleEndian) 487 { 488 enum MANTISSA_LSB = 0; 489 enum MANTISSA_MSB = 1; 490 } 491 else 492 { 493 enum MANTISSA_LSB = 1; 494 enum MANTISSA_MSB = 0; 495 }